NEFSC Trawl Bodymass Distribution

Patterns in Biomass Allocation

The following code digs into where the biomass is allocated among the 62 species sampled by the NE Groundfish Survey. Stratified abundances are used. Individual lengths and bodymasses are used to isolate what proportion of the overall abundances and biomasses reside.

Data is prepared and updated using {targets} to ensure a consistent data state and a reproducible workflow.

Target Data

Data for the report comes directly from the {targets} workflow found in _targets.R. For this markdown I’m loading in the results from the size spectrum slope analysis and the data that went into it so I can dig into any odd patterns I see. There is also regional SST based on the same trawl regions.

####  Data  ####

####  Size Spectrum Targets 

# SS and manual bins together
withr::with_dir(rprojroot::find_root('_targets.R'), 
                tar_load(size_spectrum_indices))   

# Format Columns
size_spectrum_indices <- size_spectrum_indices  %>% 
  mutate(season = fct_rev(season),
         survey_area = factor(survey_area, 
                              levels = c("GoM", "GB", "SNE", "MAB")),
         yr = as.numeric(as.character(Year)))



####  OISST Data

# OISST for all the regions
withr::with_dir(rprojroot::find_root('_targets.R'), 
                tar_load(regional_oisst))          


####  Size Data Targets

# 1. Biological data used as input
withr::with_dir(rprojroot::find_root('_targets.R'), 
                tar_load(nefsc_1g))            

# quick little format
nefsc_1g <- nefsc_1g %>% 
  mutate(fishery = case_when(
    fishery == "com" ~ "Commercially Targeted",
    fishery == "nc" ~ "Not Commercially Targeted",
    TRUE ~ "Not Labelled"))

Exploration of Abundance/Size Information

The biological data was filtered prior to doing the size spectrum analysis so that any individuals smaller than 1g were removed. Then to dig into where the biomass andd abundances were distributed I’ve made some size bins to help tease out what the community looks like.

# Cut up discrete length and weight bins
nefsc_size_bins <- nefsc_1g %>% 
  mutate(
    length_bin = case_when(
      length_cm <= 5   ~ "0 - 5cm",
      length_cm <= 10  ~ "5 - 10cm",
      length_cm <= 25  ~ "10 - 25cm",
      length_cm <= 50  ~ "25 - 50cm",
      length_cm <= 75  ~ "50 - 75cm",
      length_cm <= 100 ~ "75 - 100cm",
      length_cm >= 100 ~ ">= 100cm"),
    length_bin = factor(length_bin, levels = c(
      "0 - 5cm", "5 - 10cm", "10 - 25cm",
      "25 - 50cm", "50 - 75cm", "75 - 100cm", ">= 100cm")))


# Weight bins
nefsc_size_bins <- nefsc_size_bins %>% 
  mutate(
    weight_bin = case_when(
      ind_weight_kg <= 0.001 ~ "0 - 1g",
      ind_weight_kg <= 0.005 ~ "1 - 5g",
      ind_weight_kg <= 0.010 ~ "5 - 10g",
      ind_weight_kg <= 0.050 ~ "10 - 50g",
      ind_weight_kg <= 0.100 ~ "50 - 100g",
      ind_weight_kg <= 0.500 ~ "100 - 500g",
      ind_weight_kg <= 1.000 ~ ".5 - 1kg",
      ind_weight_kg <= 5.000 ~ "1 - 2kg",
      ind_weight_kg <= 5.000 ~ "2 - 5kg",
      ind_weight_kg <= 10.00 ~ "5 - 10kg",
      ind_weight_kg >= 10.00 ~ ">= 10kg" ),
    weight_bin = factor(weight_bin, levels = c(
      "0 - 1g", "1 - 5g", "5 - 10g", "10 - 50g",
      "50 - 100g", "100 - 500g", ".5 - 1kg",
      "1 - 2kg", "2 - 5kg", "5 - 10kg", ">= 10kg")))


# Rename the functional groups
nefsc_size_bins <- nefsc_size_bins %>% 
  mutate(
    spec_class = case_when(
      spec_class == "gf"  ~ "Groundfish",
      spec_class == "pel" ~ "Pelagic",
      spec_class == "dem" ~ "Demersal",
      spec_class == "el"  ~ "Elasmobranch",
      spec_class == "<NA>" ~ "NA"))

# Make regions go N->S
nefsc_size_bins <- nefsc_size_bins %>% 
  mutate(survey_area = factor(survey_area, levels = c("GoM", "GB", "SNE", "MAB")),
         season = factor(season, levels = c("Spring", "Fall")))

Group Summary Functions

The following function is used to process the same numbers for each combination of factor groups without me rewriting the function over and over again.

# Function to process summaries for various factor combinations
get_group_summaries <- function(...){
  
  # Do some grouping to get totals
  group_totals <- nefsc_size_bins %>% 
    group_by(...) %>% 
    summarise(total_survey_catch = sum(numlen, na.rm = T),
              total_lw_bio = sum(sum_weight_kg, na.rm = T),
              total_strat_abund = sum(strat_total_abund_s, na.rm = T),
              total_strat_lw_bio = sum(strat_total_lwbio_s, na.rm = T), 
              .groups = "drop") 
  
  # length bins
  group_lengths <- nefsc_size_bins %>% 
    group_by(..., length_bin) %>% 
    summarise(lenbin_survey_catch = sum(numlen),
              lenbin_lw_bio       = sum(sum_weight_kg),
              lenbin_strat_abund  = sum(strat_total_abund_s),
              lenbin_strat_lw_bio = sum(strat_total_lwbio_s), 
              .groups = "drop") %>% 
    left_join(group_totals) %>% 
    mutate(
      perc_total_catch  = (lenbin_survey_catch - total_survey_catch) * 100,
      perc_lw_bio       = (lenbin_lw_bio - total_lw_bio) * 100,
      perc_strat_abund  = (lenbin_strat_abund - total_strat_abund) * 100,
      perc_strat_lw_bio = (lenbin_strat_lw_bio - total_strat_lw_bio) * 100)
  
  # weight bins
  group_weights <- nefsc_size_bins %>% 
    group_by(..., weight_bin) %>% 
    summarise(wtbin_survey_catch = sum(numlen),
              wtbin_lw_bio       = sum(sum_weight_kg),
              wtbin_strat_abund  = sum(strat_total_abund_s),
              wtbin_strat_lw_bio = sum(strat_total_lwbio_s),
              .groups = "drop") %>% 
    left_join(group_totals) %>% 
    mutate(
      perc_total_catch  = (wtbin_survey_catch / total_survey_catch) * 100,
      perc_lw_bio       = (wtbin_lw_bio / total_lw_bio) * 100,
      perc_strat_abund  = (wtbin_strat_abund / total_strat_abund) * 100,
      perc_strat_lw_bio = (wtbin_strat_lw_bio / total_strat_lw_bio) * 100)

return(list("length_bins" = drop_na(group_lengths),
            "weight_bins" = drop_na(group_weights)))
}


# Process each group
regions <- get_group_summaries(Year, survey_area)
seasons <- get_group_summaries(Year, season)
fgroups <- get_group_summaries(Year, spec_class)
fishery <- get_group_summaries(Year, fishery)

Same idea here, but for plotting them based on their groupings

# 1. Abundance by length
abundance_by_length <- function(length_summary, facet_var){

  ggplot(length_summary,
       aes(Year, lenbin_strat_abund, fill = length_bin)) +
  geom_col(position = "stack", color = "gray80", size = 0.1) +
  facet_wrap(as.formula(paste("~", facet_var)), scales = "free") +
  scale_fill_gmri() +
  scale_y_continuous(labels = comma_format()) +
  labs(y = "Stratified Total Abundance",
       title = "Abundance Allocation by Length",
       x = "",
       fill = "Individual Length (cm)")



}

# test
# abundance_by_length(regions$length_bins, "survey_area")


# 2. Biomass by length
biomass_by_length <- function(length_summary, facet_var){

ggplot(length_summary,
       aes(Year, lenbin_strat_lw_bio, fill = length_bin)) +
  geom_col(position = "stack", color = "gray80", size = 0.1) +
  facet_wrap(as.formula(paste("~", facet_var)), scales = "free") +
  scale_fill_gmri() +
  scale_y_continuous(labels = comma_format()) +
  labs(y = "Stratified Total Biomass (kg)",
       title = "Biomass Allocation by Length",
       x = "",
       fill = "Individual Length (cm)")

}


# 3. Abundance by bodymass
abundance_by_bodymass <- function(bodymass_summary, facet_var){

  ggplot(bodymass_summary,
         aes(Year, wtbin_strat_abund, fill = weight_bin)) +
    geom_col(position = "stack", color = "gray80", size = 0.1) +
    facet_wrap(as.formula(paste("~", facet_var)), scales = "free") +
    scale_fill_gmri() +
    scale_y_continuous(labels = comma_format()) +
    labs(y = "Stratified Total Abundance",
         title = "Abundance Allocation by Bodymass",
         x = "",
         fill = "Individual Weight (kg)")


}





# 4. Biomass across bodymass
biomass_by_bodymass <- function(bodymass_summary, facet_var){

  ggplot(bodymass_summary,
         aes(Year, wtbin_strat_lw_bio, fill = weight_bin)) +
    geom_col(position = "stack", color = "gray80", size = 0.1) +
    facet_wrap(as.formula(paste("~", facet_var)), scales = "free") +
    scale_fill_gmri() +
    scale_y_continuous(labels = comma_format()) +
    labs(y = "Stratified Total Biomass (kg)",
         title = "Biomass Allocation by Bodymass",
         x = "",
         fill = "Individual Weight (kg)")


}

Regional Differences

These first plots get at how the different regions compare to one another. Are we seeing stiking patterns across them all simultaneously, or are there localized patterns that only occur in some of them.

# Do some grouping
regions <- get_group_summaries(Year, survey_area)

Length Bins

# Get fraction for different length/weight classes

# length bins
regional_lengths <- regions$length_bins

# Abundance by length
abundance_by_length(regional_lengths, "survey_area")

# biomass
biomass_by_length(regional_lengths, "survey_area")

Weight Bins

# weight bins
regional_weights <- regions$weight_bins

abundance_by_bodymass(regional_weights, "survey_area")

biomass_by_bodymass(regional_weights, "survey_area")

Seasonal Differences

# Do some grouping
seasons <- get_group_summaries(Year, season)

Length Bins

# Get fraction for different length/weight classes

# length bins
seasonal_lengths <- seasons$length_bins

abundance_by_length(seasonal_lengths, "season")

biomass_by_length(seasonal_lengths, "season")

Weight Bins

# weight bins
seasonal_weights <- seasons$weight_bins

abundance_by_bodymass(seasonal_weights, "season")

biomass_by_bodymass(seasonal_weights, "season")

Functional Group Differences

# Do some grouping
fgroups <- get_group_summaries(Year, spec_class)

Length Bins

# Get fraction for different length/weight classes

# length bins
fgroup_lengths <- fgroups$length_bins

abundance_by_length(fgroup_lengths, "spec_class")

biomass_by_length(fgroup_lengths, "spec_class")

Weight Bins

# weight bins
fgroup_weights <- fgroups$weight_bins

abundance_by_bodymass(fgroup_weights, "spec_class")

biomass_by_bodymass(fgroup_weights, "spec_class")

Functional Group Key

Thought it may be helpful to have a table here to check which species are currently assigned to which group.

# Display table of whichever specis fall into the categories
nefsc_size_bins %>% 
  distinct(comname, scientific_name, spec_class) %>% 
  rename(`Common Name` = comname, 
         `Scientific Name` = scientific_name,
         `Functional Group` = spec_class) %>% 
  arrange(`Functional Group`, `Common Name`) %>% 
  DT::datatable(filter = "top", rownames = FALSE)

Fishery Status Diffferences

# Do some grouping
fishery <- get_group_summaries(Year, fishery) 

Length Bins

# Get fraction for different length/weight classes

# length bins
fish_lengths <- fishery$length_bins

abundance_by_length(fish_lengths, "fishery")

# # Biomass from Survey using L-W
biomass_by_length(fish_lengths, "fishery")

Weight Bins

# weight bins
fish_weights <- fishery$weight_bins

abundance_by_bodymass(fish_weights, "fishery")

biomass_by_bodymass(fish_weights, "fishery")

Fishery Group Key

Thought it may be helpful to have a table here to check which species are currently assigned to which group.

# Display table of whichever specis fall into the categories
nefsc_size_bins %>% 
  distinct(comname, scientific_name, fishery) %>% 
  rename(`Common Name` = comname, 
         `Scientific Name` = scientific_name,
         `Fishery Status` = fishery) %>% 
  arrange(`Fishery Status`, `Common Name`) %>% 
  DT::datatable(filter = "top", rownames = FALSE)

---

Community Size Spectra Results

Community size spectrum slopes were estimated using 2 methods for comparison. An individual size distribution approach ISD developed by A. Edwards and using a more traditional method of binning biomass information into bins before fitting a slope to those bins.

The first method usezthe {sizeSpectra} package which were shown to be the most accurate when using simulated data compared to any of the binning methods.

The second method uses binned abundances, with bodymass bins of width 0.5 on a log10 scale, so 10^0 - 10^0.5 etc. These bins were then normalized by dividing the abundances by the bin width to account for the increasing bin width.

Starting data used for both methods was the same. A minimum individual biomass of 1 gram was used to avoid issue with sampling biases for smaller individuals. Area-stratified total abundances (and their corresponding length-based biomasses) were used to preserve the importance of the sampling design.

Pull Results from Both Methods

# Pull the group ID for the slopes grouped on year, season, and region
warmem_group_slopes <- size_spectrum_indices %>% 
  filter(`group ID`== "single years * season * region")


# Or just regions and years
year_region_slopes <- size_spectrum_indices %>% 
  filter(`group ID` == "single years * region")

Biomass Data to Match Size Spectrum Analysis

Region & Season

The first lens to look at is the seasonal variation across all the different areas. This was the typical grouping that we have been focusing on.

# Do some grouping
warmem_totals <- get_group_summaries(Year, survey_area, season)


ggplot(warmem_totals$weight_bins, 
       aes(Year, wtbin_strat_lw_bio, fill = weight_bin)) +
  geom_col(position = "stack", color = "gray80", size = 0.1) +
  facet_grid(survey_area~season, scales = "free") +
  scale_fill_gmri() +
  scale_y_continuous(labels = comma_format()) +
  labs(y = "Stratified Total Biomass (kg)",  
       title = "Biomass Allocation by Weight",
       x = "", 
       fill = "Individual Weight (kg)")

Region & Functional Group

The second lens that I feel is important is the functional groups. This is the not the typical grouping that we have been focusing on, but clears the air on where the biomass increase is coming from.

# Do some grouping
fgroup_area <- get_group_summaries(Year, survey_area, spec_class)

ggplot(fgroup_area$weight_bins, 
       aes(Year, wtbin_strat_lw_bio, fill = weight_bin)) +
  geom_col(position = "stack", color = "gray80", size = 0.1) +
  facet_grid(survey_area~spec_class, scales = "free") +
  scale_fill_gmri() +
  scale_y_continuous(labels = comma_format()) +
  labs(y = "Stratified Total Biomass (kg)",  
       title = "Biomass Allocation by Weight",
       x = "", 
       fill = "Individual Weight (kg)")

Edwards Methodology

The Edwards methodology differs from the other methods for estimating size spectra by avoiding the subjective decisions around how to bin data prior to fitting the log-linear regression.

Instead, Edwards uses the individual size distributions (how many individuals of any given length/weight). Abundances are totaled into discrete size bins based on expected biomass at length and length + 1, and individuals that fall within those bins are totaled to get abundances across a continuous distribution of individual bodymass.

The next difference is that the individual body size data is presumed to follow a bounded power-law distribution, with a minimum and maximum body size. Using the individual size data, maximum likelihood estimation is used to solve for the parameter (b) that minimizes the negative log-likelihood.

Year and Region Only

# Just years and regions
(ss_patterns <- year_region_slopes %>% 
  ggplot(aes(yr, b, color = survey_area)) +
  geom_line(aes(group = survey_area), linetype = 3) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "lm", formula = y ~ x, 
              alpha = 0.2) +
  facet_grid(survey_area~.) +
  scale_color_gmri() +
  labs(x = "", 
       y = "Size Spectrum Slope (b)", 
       color = "",
       title = "Results from Area-Stratified Abundance - Edwards Method"))

Year, Region, and Season

# Plot the sizeSpectra slopes - year*region*season
(ss_patterns <- warmem_group_slopes %>% 
  ggplot(aes(yr, b, color = survey_area)) +
  geom_line(aes(group = survey_area), linetype = 3) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "lm", formula = y ~ x, 
              alpha = 0.2) +
  facet_grid(survey_area~season) +
  scale_color_gmri() +
  labs(x = "", 
       y = "Size Spectrum Slope (b)", 
       color = "",
       title = "Results from Area-Stratified Abundance - Edwards Method") 
)

Odd Years

In the previous figure, for Georges Bank there are a number of years when b is close to -1.2 that do not follow the rest of the trend. These occur on 1974 and 1980. This section exists to see why that is:

# Find the years
odd_years <- warmem_group_slopes %>% 
  filter(survey_area == "GB", 
         b < -1.2) %>% 
  split(.$Year)

# Find the data
odd_year_dat <- nefsc_size_bins %>% 
  filter(survey_area == "GB",
         Year %in% c(1974, 1980),
         season == "Fall") %>% 
  split(.$Year)

Individual Size Distribution

The individual size distribution plots look like alot but are not too complicated with a little explanation.

The blue bars are the width of what a given species biomass can be based on. What the species is, and its L-W relationship. Since fishes are measured to the nearest cm the leftbound is the biomass for that length, and the right bound is how much mass it could also have before jumping into the next 1cm increment.

The height of the blue bars indicate how many individuals fall within that size range. Eventually building up the range of possible sizes caught, and how many there are/were.

The red line is a bounded power law fit to those blue lines. The power law relationship between abundance and size is the foundation for all the size spectra methods, in this case individual sizes are used rather than the more typical larger bins.

Data Prep

# ISD plot preparation
# counts how many of each size there are with overlapping sizes counting
# towards the discrete size bins when possible
isd_prepped <- map(odd_year_dat, ~ isd_plot_prep(as.data.frame(.x),
                                                stratified_abundance = T,
                                                min_weight_g = 1))

Plotting Function

# Function for ISD plots 
# should work on whatever the grouping variable is as long as the lists match
#' @title Plot Individual Size Distribution Curves
#'
#' @param isd_data_prepped Data from isd_plot_prep
#' @param mle_results Corresponding group results from group_mle_calc
#' @param abundance_used Label to add for context of what abundance source was used
#' @param plot_rects Flag to plot bodymass rectangles
#' @param show_pl_fit Flag to include powerlaw fit
#' @param xlim_global global xlimits so that other plots can be held side-by-side
#' @param group_name name of the group for the data to use as label
#'
#' @return
#' @export
#'
#' @examples
ggplot_isd <- function(isd_data_prepped, 
                       mle_results, 
                       abundance_vals = "observed",
                       plot_rects = TRUE,
                       show_pl_fit = TRUE,
                       xlim_global = NULL,
                       group_name = NULL){
  
  # Columns and labels
  abundance_lab <- c("observed" = "Survey Abundance", "stratified" = "Stratified Abundance")
  abundance_label <- abundance_lab[[abundance_vals]]
  
  # Power law parameters and summary details for the group of data:
  b.MLE           <- mle_results$b
  total_abundance <- mle_results$n
  b.confMin       <- mle_results$confMin
  b.confMax       <- mle_results$confMax
  
  
  # Create range of x values from the group to get power law predictions
  # PLB = bounded power-law
  # min and max weights for predictions
  xmin  <- mle_results$xmin
  xmax  <- mle_results$xmax
  
  # Create x values (individual bodymass) to predict across
  # break up the Xlim into pieces between min and max
  x.PLB <- seq(from = xmin, 
               to   = xmax, 
               length.out = 2000)   
  
  # get the length of that vector
  x.PLB.length <- length(x.PLB)  
  
  # remove last entry, add an entry .9999 of the way there, and cap it with the last entry wtf
  x.PLB <- c(x.PLB[-x.PLB.length], 0.99999 * x.PLB[x.PLB.length], x.PLB[x.PLB.length])
  
  
  # Y values for plot limits/bounds/predictions from bounded power law pdf
  y.PLB         <- (1 - sizeSpectra::pPLB(x = x.PLB, b = b.MLE, xmin = min(x.PLB), xmax = max(x.PLB))) * total_abundance
  y.PLB.confMin <- (1 - sizeSpectra::pPLB(x = x.PLB, b = b.confMin, xmin = min(x.PLB), xmax = max(x.PLB))) * total_abundance
  y.PLB.confMax <- (1 - sizeSpectra::pPLB(x = x.PLB, b = b.confMax, xmin = min(x.PLB), xmax = max(x.PLB))) * total_abundance
  
  
  # Put it in a df to make it easier
  PLB_df <- data.frame(
    x.PLB   = x.PLB,
    y.PLB   = y.PLB,
    confMin = y.PLB.confMin,
    confMax = y.PLB.confMax)
  
  
  # Coefficient Labels
  group_name <- str_replace_all(group_name, "_", " ")
  the_slope  <- as.character(round(mle_results$b, 3))
  
  
  ####  First Plot  
  p1 <- ggplot()
  
  # Toggle Rectangles
  if(plot_rects == TRUE) {
    p1 <- p1 + geom_rect(data = isd_data_prepped, 
                         aes(xmin = wmin_g, 
                             xmax = wmax_g, 
                             ymin = lowCount, 
                             ymax = highCount),
                         color = "gray70",
                         fill = "transparent")}
  
  # Add segments for bin widths
  p1 <- p1 + geom_segment(data = isd_data_prepped,
                          aes(x = wmin_g, 
                              xend = wmax_g, 
                              y = countGTEwmin, 
                              yend = countGTEwmin),
                          color = "blue")
  
  
  # Set limits to global limits if desired
  if(is.null(xlim_global) == FALSE){
    p1 <- p1 + scale_x_log10(limits = xlim_global, 
                             labels = trans_format("log10", math_format(10^.x)))
  } else{
    p1 <- p1 + scale_x_log10(labels = trans_format("log10", math_format(10^.x)))
  }
  
  
  # Toggle to turn on the fit lines or not
  if(show_pl_fit == TRUE) {
    p1 <- p1 +
      geom_line(data = PLB_df, aes(x.PLB, y.PLB), color = "darkred") +
      geom_line(data = PLB_df, aes(x.PLB, confMin), color = "darkred", linetype = 2) +
      geom_line(data = PLB_df, aes(x.PLB, confMax), color = "darkred", linetype = 2) }
  
  # Finish off labels
  p1 <- p1 + labs(x = "Individual Bodymass (g)",
                  y =  "Fishes with BodyMass \u2265 x",
                  title = group_name, 
                  tag = str_c("b = ", the_slope),
                  caption = str_c(abundance_label, " used for size spectrum slope estimation.")) +
    theme(plot.tag.position = "topright")
  
  

  ####  Second Plot, log transformed y axis
  
  # Start plot
  p2 <- p1
  
  # Change the y-axis, log transformation
  p2 <- p2 + scale_y_log10(labels = trans_format("log10", math_format(10^.x))) +
    labs(x = "Individual Bodymass (g)",
         y =  "Fishes with BodyMass \u2265 x",
         title = group_name, 
         tag = str_c("b = ", the_slope),
         caption = str_c(abundance_label, " used for size spectrum slope estimation.")) +
    theme(plot.tag.position = "topright")
  
  
  # Build the stacked plot
  p3 <- p1 / p2
  
  # Put figures in list
  plot_list <- list(obs_y   = p1,
                    log10_y = p2,
                    stacked = p3)
  
  return(plot_list)
}

Individual Size Distributions

Like any of these methods sometimes the assumed relationship does not perfectly fit the data. In the case of these two years there appears to be in the 10g-1kg range.

isd_1974 <- ggplot_isd(isd_data_prepped = isd_prepped$`1974`, 
                       mle_results = odd_years$`1974`, 
                       abundance_vals = "stratified", 
                       plot_rects = FALSE, 
                       show_pl_fit = TRUE, 
                       group_name = "1974, Fall, GB")$obs_y

isd_1980 <- ggplot_isd(isd_data_prepped = isd_prepped$`1980`, 
                       mle_results = odd_years$`1974`, 
                       abundance_vals = "stratified", 
                       plot_rects = FALSE, 
                       show_pl_fit = TRUE, 
                       group_name = "1980, Fall, GB")$obs_y

isd_1974 / isd_1980

Binned Biomass

If I plot the same data using the binning approach you still see that the bins that cover 10^1 through 10^3 grams are higher than the fit line (analogous to red line for ISD). But, because the size distribution is no longer continuous, you can potentially run into issues of bins not accurately representing the data that goes into them.

# Log10 binning, 0.5 increments on log10 scale
# assign bins
l10_assigned <- map(odd_year_dat, ~ assign_log10_bins(.x))

# plot bin structure for survey abundance and stratified abundance
l10_1974 <- plot_log10_ss(l10_assigned = l10_assigned$`1974`) 
l10_1980 <- plot_log10_ss(l10_assigned = l10_assigned$`1980`) 


# Grab just the stratified abundance one
l10_1974[[2]] + labs(title = "1974, Fall, GB") |
l10_1980[[2]] + labs(title = "1980, Fall, GB")

Log10 Binning

Slope

# plot trends of log10 slopes
(l10_patterns <- warmem_group_slopes %>% 
  ggplot(aes(yr, l10_slope_strat, color = survey_area)) +
  geom_line(aes(group = survey_area), linetype = 3) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "lm", formula = y ~ x, 
              alpha = 0.2) +
  facet_grid(survey_area~season) +
  scale_color_gmri() +
  labs(x = "", 
       y = "Size Spectrum Slope (b)", 
       color = "",
       title = "Results from Area-Stratified Abundance - log10 bins"))

Intercept

# plot trends of log10 slopes
(int_patterns <- warmem_group_slopes %>% 
  ggplot(aes(yr, l10_int_strat, color = survey_area)) +
  geom_line(aes(group = survey_area), linetype = 3) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "lm", formula = y ~ x, 
              alpha = 0.2) +
  facet_grid(survey_area~season) +
  scale_color_gmri() +
  labs(x = "", 
       y = "Size Spectrum Intercept", 
       color = "",
       title = "Results from Area-Stratified Abundance - log10 bins")
)

Fit Statistics (adjusted r-squared)

For each of the slope/intercepts derived using a linear model and binned data I also pulled the adjusted r-square to get a sense of whether or not certain groups had poor fits that should be investigated.

# Reshaping
# Goal: Row for Years, columns for each different group
l10_fit_dat <- warmem_group_slopes %>% 
  select(Year, season, survey_area, l10_slope_strat, l10_int_strat, l10_rsqr_strat) %>% 
  mutate(yr = as.numeric(Year))


l10_fit_dat %>% 
  ggplot(aes(yr, l10_rsqr_strat, color = survey_area)) +
  geom_rect(xmin = -Inf, xmax = Inf, ymin = 0, ymax = 0.2, 
            fill = "gray80", color = "transparent") +
  geom_hline(yintercept = 0.2, linetype = 2, size = 0.5, color = "gray50") +
  geom_line(aes(group = survey_area), linetype = 3) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "lm", formula = y ~ x, 
              alpha = 0.2) +
  facet_grid(survey_area~season) +
  scale_color_gmri() +
  scale_y_continuous(limits = c(0,1), expand = c(0,0)) +
  labs(x = "", 
       y = "Linear Model R-Square", 
       color = "",
       title = "Results from Area-Stratified Abundance - log10 bins")

Sea Surface Temperature

Long-Term Patterns

# plot trends of log10 slopes
(sst_patterns <- regional_oisst %>% 
  ggplot(aes(yr, sst_anom, color = survey_area)) +
  geom_line(aes(group = survey_area), linetype = 3) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "lm", formula = y ~ x, 
              alpha = 0.2) +
  facet_grid(survey_area~.) +
  scale_color_gmri() +
  labs(x = "", 
       y = "Mean Temperature Anomaly", 
       color = "",
       title = "")
)

 

---

A work by Adam A. Kemberling

Akemberling@gmri.org